Forming doublons by a quantum quench 
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Repulsive interactions between particles on a lattice may lead to bound states, so called doublons. 
Such states may be created by dynamically tuning the interaction strength, e.g. using a Feshbach 
resonance, from attraction to repulsion. We study the doublon production efficiency as a function 
of the tuning rate at which the on-site interaction is varied. An expectation based on the Landau- 
Zener law suggests that exponentially few doublons are created in the adiabatic limit. Contrary to 
such an expectation, we found that the number of produced doublons scales as a power law of the 
tuning rate with the exponent dependent on the dimensionality of the lattice. The physical reason 
for this anomaly is the effective decoupling of doublons from the two-particle continuum for center 
of mass momenta close to the corners of the Brillouin zone. The study of doublon production may 
be a sensitive tool to extract detailed information about the band structure. 



Introduction - Quantum evolution of a system driven 
far from equilibrium by time-dependent perturbations 
became a focus of increased attention with the advent 
of cold atomic gases. These systems provide a unique 
advantage due to their high degree of tunability vis-a-vis 
various system parameters including dimensionality, den- 
sity, inter-particle interaction strengths and disorder [T]. 
In particular, such systems provide remarkably accurate 
realizations of both Bose and Fermi Hubbard models, 
in which the rich equilibrium phase diagrams [5] or far- 
from-equilibrium dynamics may be studied in a highly 
controllable manner. 

For strong attraction, spin— 1/2 fermions preferentially 
pair on lattice sites forming tightly bound composite 
bosons. A slow quench of the interaction strength takes 
the system through the BEC-BCS crossover into a regime 
of strong repulsion, resulting in a ground state which is 
a Fermi liquid or Mott insulator, depending on the fill- 
ing fraction However, if such a quench is sufRciently 
fast, a fraction of initial pairs may not have time to dis- 
sociate and thus transform into repulsively bound states 
— doublons (6|. Although being high-energy states, dou- 
blons have a long lifetime due to energy conservation, 
which requires a coherent multi-particle excitation to in- 
duce doublon decay. The lifetime of residual doublons 
was investigated for both bosonic and fermionic coun- 
terparts [U [5] . Here we focus on the creation efhciency 
of such repulsively-bounded pairs upon a quench 

of the on-site interaction strength. The latter may be 
achieved, using e.g. a Feshbach resonance [7], and is 
characterized by a rate A dictating the speed of the ramp 
for the on-site particle interaction U = Xt. After the 
quench a fraction of remaining double occupancies, or 
production efficiency of doublons, depends sensitively on 
the ramp rate: for fast ramps most doublons survive, 
while in the adiabatic limit the majority softly dissociate 
into free particles states forming a Bloch band. 

Our main results for the doublon production efficiency 
are as follows. In the limit of small initial concentration 
of attractively bound pairs the problem may be treated 
in the spirit of a Landau-Zener crossing [HHTU] of a time- 




FIG. 1. (Color online) Two-particle spectrum and doublon 
bound state energy as a function of the center of mass quasi- 
momentum K for negative on-site interaction U. Due to 
quasi-translational invariance, doublons evolve along curves 
of constant K and survive the ramp with probability V{K). 
Near the corners of the Brillouin zone, \K\a = tt, the two- 
particle bandwidth and doublon-band coupling tend to zero, 
rendering those states most likely to surive the ramp. 

dependent level and a time-independent band. Such a 
variant of the Landau-Zener problem is known as the 
Demkov-Osherov model [101 E] i which admits an exact 
solution. In terms of the rate A and nearest-neighbor hop- 
ping Ji along the primitive vectors of a 13— dimensional 
lattice, we found for the probability for the two parti- 
cles to remain on the same site after a quench {e.g. for 
t +00) 

V = Y[loi8njf/\)e-'-'?/\ 

4=1 



where J — {Yli Ji)" is the geometric mean hopping, Iq 
is a modified Bessel function of the first kind and the 
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asymptotic limits are provided in the second line. Re- 
markably, the doublon survival probability is not ex- 
ponentially suppressed in the adiabatic limit, like V ^ 
g-87rj; /A^ but rather scales as a power law V ^ 
(A/J^)^/^. This implies that doublons have a signifi- 
cantly better chance of surviving the quench than might 
be naively expected. As explained below, the power- 
law survival enhancement comes from the presence of 
doublons with center of mass quasi-momenta near the 
Brillioun zone corners, see Fig. [l] We also considered 
the case of a finite concentration of initially attractive 
bounded pairs, taking into consideration the bosonic or 
fermionic nature of dissociated particles. We show that 
particle statistics renormalizes the production efficiency 
in an intuitive way: for fermions, Pauli blocking of dis- 
sociation channels leads to doublon production enhance- 
ment, while for bosons stimulated dissociation suppresses 
doublon production. 

Dilute limit - In the presence of a sufficiently deep 
optical lattice one may neglect higher Bloch bands and 
focus entirely on the lowest. This approximation leads 
to the Hubbard model characterized by hopping Ji along 
the primitive vectors and on-site interaction energy U. 
Below we express the corresponding Hamiltonian assum- 
ing particles are spin— ^ fermions. The spinless bosonic 
counterpart is obtained by removing spin indices and cor- 
responding summations. 

H = - + Uit)Y^ n,^n,^i - ^ ^ (2) 

{i,j)cr i i 

where (i, j) restricts summation to nearest neighbor sites, 
a — ±1 denotes the spin, (c) are fermionic creation 
(annhilation) operators and nja = CjaCja is the number 
operator. Below we focus on a rectangular lattice with 
lattice constants and hopping (i = 1, 2, 3). 

In the case of two particles one may write a wavefunc- 
tion of their coordinates \l/(xi, X2) in terms of central and 
relative coordinates as \l/(xi,X2) = e'-'^^'0K(r), where 
R = (xi + X2)/2, r = Xi — X2 and K is the center of 
mass quasi-momentum of the particle pair. In the case 
of bosons, one requires the relative wavefunction ^p-K^r) 
to be an even function of the relative separation r. The 
same is also true of spin— ^ fermions occupying the spin- 
singlet state (on the triplet manifold ipK{r) is odd and 
thus on-site interaction is irrelevant). Consequently, the 
two-particle results derived below extend to either species 
of constituent particles. We shall find, unsurprisingly, 
that particle statistics do play a role in the many-body 
situation considered later. 

Due to translation invariance of lattice systems, we 
work in a sector with fixed K. The action of Hamil- 
tonian ^ on the two-particle wavefunction e**^^'0K(r) 
produces an effective tight-binding model in the relative 
coordinate r 

{-3-Ar + EK + USrfi) V'K (r) - E^bK (r) , (3) 

where we introduced an effective hopping with com- 
ponents Ji = 2JiC0S^^, the discrete Laplacian 
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FIG. 2. (Color online) Top panel: average two-particle energy 
of a dissociated doublon as a function of the center of mass 
quasi-momentum K. Each curve is plotted with a different 
value of the parameter A/8J^ given by (from top to bottom): 
6, 3, 1.5, 1, 0.5, 0.2, 0.01. The dashed curve is the lower 
edge of the two-particle continuum. Bottom panel: average 
dissociation energy integrated over K as a function of the 
rate. Dashed lines correspond to the asymptotic limits given 
by Eq. Q. 

[ArV'K(r)], = 'i/'K(r + a,) -I- ?/'k(i" - a^) - 2ipic{r) and the 
center of mass kinetic energy Ek = 4Ji(l — cos^^i^). 

An interesting feature of Eq. ([s]) is that it supports a 
bound state solution at both negative and positive val- 
ues of the on-site interaction U. This can be seen most 
easily introducing the Fourier transformation 'i/'K(r) = 
giqr^j^^q-j -j^^Q gq_ which gives rise to an equa- 
tion for the bound state energy E^ 

1=/^ ^ . (4) 

J (27r)^i?d-e(K,q) ^ ' 

Here e(K, q) = ^ - 4Ji [l — cos^^i^cosfq^ai)] is the two- 
particle energy spectrum, while Eq. m may be viewed 
as a condition of the dressed Green function having a 
pole on the real axis at E — E^ ^ [H]. In one spa- 
tial dimension {e.g., 3^ — 3y = 0), the particular 
form of the two-particle spectrum allows to integrate 
Eq. Q explicitly, giving the bound state energy Ed{K) = 
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4J+sgn(f7)yC/2 + le.Pcos'^^, seen in Fig.jl] The exis- 
tence of a bound state for any value of C/ 7^ extends also 
to Z? = 2, where the density of states Af{e) near the band 
bottom is constant JV{e) — > Aq. Similar to the Cooper 
problem, the bound state energy E^i —Je~^^^°^ is 
thus exponentially shallow as ?7 — 0. In I? = 3 the den- 
sity of states vanishes at the band bottom, resulting in 
the convergence of the momentum integral at i^d = and 



thus leads to a critical value UcyH (K) — j 1^ ^^y^ 

such that for \U\ < L'crit(K) no bound state exists for a 
given K. 

Another intriguing property of Eq. (|3| is the depen- 
dence of the effective hopping on the center of mass 
momentum: J.^ ~ 2 J^cos (Kiai/2). Consequently, the 
hopping along a particular direction is identically zero 
on the corresponding Brillouin zone surface defined by 
|Ki| = Tr/ui. This is the result of destructive interference 
between two amplitudes that both increase the relative 
coordinate r along the specified direction, with one shift- 
ing the central coordinate by +ai/2, the other by —ai/2. 
The presence of zero hopping is also manifest in the two- 
particle band spectrum (see Fig. [I] for 13 — 1) which 
becomes pinched at the BZ boundaries, leading to a sys- 
tem of decoupled (and hence degenerate) states. Such 
momentum dependent hopping was discussed for D = 1 
in Refs. [HI [E] in the context of Feshbach molecules 
coupled to a band continuum. 

We proceed by transforming Eq. ([s]) into a form al- 
lowing the results of the Demkov-Osherov (DO) model 
to be conveniently applied. This is achieved by intro- 
ducing a unitary transformation '0k (r) — ■(/'K(r)5r,o + 
X]q^r,q'0K,q which diagonalizcs Eq. (j3| in the subspace 
of states having r ^ (i.e., the unitary transformation 
satisfies boundary condition Uq ^ = U^ q = 0), 



Xt 



K,q 




V^k(O) 

V'K.q 



E 



V^k(O) 



(5) 



where Jk,^ = —2^ - J^cos (Kia,i/2) Wa.^q and eK,q is the 
corresponding set of band energies. Since we deal with a 
single site defect in an infinite system, the spectrum £K,q 
coincides with the two-particle spectrum e(K, q) given 
above. Indeed, the most important aspects of the site 
defect, its energy U = Xt and coupling to band states, 
has been included explicitly. We note that, owing to 
the point boundary condition, the corresponding eigen- 
functions V'K.q are not simply plane waves except for the 
special case D — 1 (see below). 

The Hamiltonian ([s]), being now of the DO form, al- 
lows to read off the exact survival probability of the r = 

state: r{K) = J]q e'^'^l"''^"!'/^ [lOlin]- Without a pre- 
cise knowledge of the unitary matrix U, such a result is 
seemingly quite useless. Remarkably, however, one can 
now determine the survival probability explicitly without 
knowing the precise form of the unitary transformation. 
Raising the product over states to summation in the ex- 
ponent, one may show using only the unitarity property 



U''U = i that 
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Within the above LZ framework the initial {t = ~oo) 
on-site energy is taken to be U/J — —00. This implies 
that initially doublons, having infinite effective mass, are 
in site eigenstates. The initial state is thus one whose 
momenta K components are uniformly distributed over 
the entire Brillouin zone [Hj. As a result, the doublon 
survival probability is given by the (evenly weighted) av- 
erage over the BZ, namely V = J^k'^O'^)- Summing 
Eq. ([6| over K gives the main result Eq. ([T]) [15 . 

After the ramp there exists, in addition to a fraction 
of surviving doublons, a fraction of dissociated, or non- 
surviving, particles occupying band states. One may ask 
what is the non-equilibrium distribution of such parti- 
cles. To this end we require the probability of doublon 
dissociation into state q: 7'q(K) = (1 — ^q) n£k<e 
where Zq = e-^-^\JK.cW>^ ^\n\. This distribution al- 
lows to study, e.g. the average energy of non-surviving 
particles in the band. For fixed K the latter is writ- 
Eq7'q(K)e(K,q)/(l-7'(K)). In 



ten as i?avg(K) 
the adiabatic limit, A 



0, 



is exponentially small. 



and as a result of the energy summation above, we 
find 7'q(K) — > Sqo, where e(K,0) is the lower edge of 
the two-particle continuum. Fig. [l] In the same limit 
V{K)^0 so that £;avg(K) 4J,(l-cosJ^) = 

(i.e., in the adiabatic limit the two-particle energy is 
simply the center of mass kinetic energy of the disso- 
ciated doublon). In the diabatic limit, A — > 00, one 
may use 7^q(K) = 27r|JK,qP/A along with H^.^e,. = 
[Wdiag{eK,q}^^]a. a ~ '^a^.a^ X^fe to calculate the 
average dissociation energy, i?avg(K) = which 
is nothing but the band center energy. To obtain results 
between the two asymptotic limits above, one must have 
access to the explicit form of the unitary matrix U^.^q. 

For simplicity of calculation, we solved for the unitary 
transformation in 13 = 1, where the point boundary con- 
dition coincides with the surface boundary condition fa- 
miliar from a particle in a box. Hence, Ur^q = 'ij~^s'm{rq) 
where q may be identified with the relative momentum 
of band particles, since the eigenfunctions are odd com- 
binations of plane waves. Taking the continuum limit 
N ^ 00, one finds the probability of dissociating into 
the range (g, q + dq), 



_ 16J^ 
A 

X Exp 



Ka 



16J2 



sin^(ga) 



,Ka 



— — cos — — \qa — smqacosgal 
A 2 



•(7) 



One notices that 'Pq{K)\q=o^±T^ — 0, implying that dou- 
blons are effectively decoupled from states with these 
relative momenta. This is unsurprising considering the 
relative group velocity Vq — dqe{K,q) vanishes at these 
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points, i.e., if the relative velocity is zero there can be 
no dissociation. Using Eq. ([t]) one may now calculate 
Eavg{K) for various A, shown in Fig. [2] 

To monitor ramp-induced heating, one may study the 
total dissociation energy E'avg — X^k ^avg(K). For arbi- 
trary D the asymptotics of this quantity lie in the range 
dictated by the A-independent limits of ii'avg(K) given 
above. For Z? = 1 we determine the leading dependence 
on A in the vicinity of these bounds shown in Fig. [2j 



4J- 




A>8J2 
A<8J2, 



(8) 



with numerical factors C — 0.966 and C — 0.448. Equa- 
tion Q implies that the system is heated by an amount 
of order J regardless of how slow the ramp is. 

Finite doublon concentration - To go beyond the above 
two-particle physics we study the secondary quantized 
Hamiltonian H = Hq + iJjnt describing the dynamics 
of doublons coupled to particles occupying band states 
(for simplicity we hereafter set the band center to be the 
zero of energy). Below we express equations for spin— ^ 
fermions, but provide results also for spinless bosons. 

Ho = U{t) ^ dl^dK - ^ 2JiC0s{qiai)cl^Cqa, (9) 

K qiCTji 

-ffint = ^^K.qd^Cf +q,CTCK_q _^ -t-h.C. (10) 

K,q,CT 

The interaction Hamiltonian Hi^t describes the dissoci- 
ation and association of (bosonic) doublons with total 
momentum K onto two particles with momenta y ± q. 
Notice that in the form Eq. (|9| , the doublon band is com- 
pletely flat E(K) — U. To obtain the renormalized dou- 
blon energy E^, one integrates out fermionic degrees of 
freedom and expands the resulting action to second order 



in c^K operators, dj^{e) 



Uk., 



dK(e) 



^ Z^q £-e(K,q) 

d^(e)2?^^(K, e)dK(£)- The doublon bound state energy 
is defined by the pole of the corresponding propaga- 
tor, I?^^(K,_Ed) = 0, which also shares the solution to 
Eq. (|4|. 

Following Ref. [T5], we derive a set of coupled quan- 
tum kinetic equations governing the distribution func- 

r{t) for 

(11) 



tions rid = nd(K,t) for doublons and n^^ 
free particles, 

nd = -7r^|JK,qp5(Ai-e(K,q)) 

q,(T 



X 



{nd (l 



-2+q,cr 2--q,-cr 



"f+q,<T^if-q 



71, 



J, = -271^ I JK^q_K/2)l''5(At - £(K, q - K/2)) 



K 



q,-<T 



- rtd (1 - "qtr - "-K-q-o-)} ■ (12) 



The factor in curly brackets in Eqs. (Ill, 



(12 



and "out" terms: 



-2-H-qo^ "2 — q. — 



(l+^d) 



rid (l - ?iK+q^) 



2 4' 



The latter piece 

corresponds to the "out" process, where 1 — 'T-K-tqio- 
reflects the fact that doublon decay is prohibited if either 
free particle state is occupied, i.e. ny^j^^j^^ = 1. For 
bosons, the minus signs in front of changed to 

plus, replacing Pauli blocking with stimulated emission. 
The "in" process may be understood on the same 
grounds. 

In the single doublon, vanishing concentration, limit 
('^qcr = 0), Eq. (11) indeed posseses the LZ solution: 

nd(K,-Hoo)/nd(K,-oo) = e-^^^-* l^'^-'l'/^, Eq. 
The aim of the above outlined kinetic scheme is thus 
to obtain leading order corrections to the efficiency 
in the presence of finite doublon concentration. In 
the rapid ramp limit, ^ A, we solve the cou- 
pled kinetic equations by first substituting the time- 
independent initial doublon distribution rid = Ua into 



Eq. ( 12 ) and solve for n„a neglecting non-linear con- 



tributions, n^„(t) = 27rn^'U ^ X^k l'^K,q-K/2p©(At - 
£(K, q— K/2)). This approximation amounts to using 

,W _ OY- .(/) , V r)^-^) 



the conservation law 2 — 2 ; 



where 



,(/) 



is calculated 



- n«(l-2^A-iEqkK,q| 
in the leading (zero density) order. Substituting the 
above approximation for n^^ into Eq. (11) using the 
D = 1 unitary matrix U above gives the hnal fraction 
of surviving doublons 



V = 1 



8n.P 
A 



4 V 4 



8nJ' 



(13) 



valid to the leading order in <C 1/2, where ± refers to 
fermionic (bosonic) particles and v is the filling fraction 

related to the initial doublon distribution as n'^^^K) = v. 
That the leading correction in v occurs in second order 
follows from the fact that one order is required to first 
populate the free-particle levels and another is needed to 
act back on the doublon decay rate. 

The more interesting (power-law) part of the survival 
probability occurs in the adiabatic regime. As explained 
in Ref. 16J, this regime is beyond the scope of the kinetic 
equation approach. In brief, this is because Eqs. (11), 



( 12 ) do not predict evolution from ground state to ground 
state, namely that all doublons must be exhausted in 
the formation of a Fermi sea as [/ — > +oo and A — 0. 
To properly take this fact into account, we propose the 
following combined strategy. In the extreme adiabatic 
limit, only those doublons with momenta components 
sufficiently close to |Ki| — n/ai have an appreciable 
chance of surviving, i.e., an effective coupling small com- 
pared to -s/A. This fact allows to treat their dynamics in 
the diabatic sense, subject to the kinetic theory outlined 
above. All other doublons with momentum away from 
the BZ boundary dissociate into states near q = Q due to 
the exponential suppression of their survival probability. 
From the perspective of doublons near the BZ corners, 
they observe an essentially uniform distribution of occu- 
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pied band states Uqa = i^- According to Eq. (11), this 



modifies the effective rate as A — >■ A/(1=f2z^). As a result, 
the power-law of Eq. ([T]) becomes 



V = 



A 



(1 =F 2v) 167r2 J2 



D/2 



(14) 



where =F refers to fermionic (bosonic) constituent parti- 
cles. From Eq. (14 1, we see that the effect of stimulated 



emission for bosons decreases the total survival probabil- 
ity, making the process appear effectively slower, while 
the effect of Pauli blocking in the case fermions increases 
the survival probability because there is less phase space 
for doublon dissociation. 

Conclusions - The presence of a periodic potential 
leads to the formation of interband gaps in which a high 
energy bound state (doublon) may exist. The formation 
of such doublons may be achieved employing the con- 



cept of Feshbach resonance in which the inter-particle 
coupling strength is changed by applying a time-varying 
magnetic field. We showed that doublon creation is not 
exponentially suppressed in the adiabatic limit but rather 
scales as a power law V ^ (A/J^)^/^, contrary to con- 
ventional LZ wisdom. Beyond the two-particle picture, 
we derived finite concentration corrections to the sur- 
vival probability depending on the bosonic or fermionic 
nature of the constituent particles. Going beyond the 
dilute limit towards filling factors of order 1/2 requires 
a more careful treatment of the Hubbard Hamiltonian 
Eq. ([2]) than provided here. In particular, the bosonic 
nature of the doublon fields c^k, rf^ in Eq. ^ can only 
be justified in the dilute limit where doublons are sparse, 
their mutual scattering is irrelevant and hence their in- 
ternal structure (composite nature) is of no importance. 
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